1 

TITLE OF THE INVENTION 

Method of Extracting Physical Model Parameter and Storage Medium Therefor, 
and Method of Manufacturing Non-linear Element 

5 BACKGROUND OF THE INVENTION 
Field of the Invention 

The present invention relates to physical properties that provide characteristic 
quantity sets g, (i = 1, 2, m) each consisting at least one characteristic quantity, 
corresponding to a plurality of extrinsic factor sets Vj, respectively, each consisting of at 
10 least one extrinsic factor, more particularly to a technique for extracting a parameter set P 
in a physical model for expressing a calculated value of each characteristic quantity set g s 
(s is any one value of i) as a function f (v s , P) of the corresponding extrinsic factor v s and 
the parameter set P consisting of a plurality of parameters. 

One example is a technique for extracting model parameters for use in a circuit 
15 simulation in LSI (Large Scale Integrated circuit) designs. 

Description of the Background Art 

A process of manufacturing an LSI is broadly divided into a design step for 
designing circuits thereof and a semiconductor processing step for implementing the 
20 circuits as a semiconductor device form according to information obtained in the design 
step. In the design step, a circuit simulation is performed to predict the functions of an 
LSI circuit as a semiconductor device (hereinafter referred to as "LSI devices"). 

We can look at the circuit simulation from two important viewpoints: 
formulation of circuit equations and device modeling. In the device modeling, electric 
25 properties of a non-linear device such as a transistor are simulated by using an analytic 
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formula obtained by modeling the device. The analytic formula includes physical or 
semi-empirically-determined parameters . 

To perform the circuit simulation with high precision, the parameters for use in 
the device modeling need to be determined properly. As an index to this determination, 
5 errors between measured characteristics of the LSI device and calculated values on the 
basis of an analytical model are usually adopted. Alternatively, instead of the measured 
characteristics of the LSI device, results of device simulation for simulating phenomena 
which occur in a device such as a transistor can be used. 

Discussion will be now on a background-art method of extracting physical 

10 model parameters. The physical model is set to obtain the physical properties of a 
non-linear element, where a non-linear relation is established between extrinsic factor sets 
Vj (i = 1, 2, 3, m) each consisting of at least one extrinsic factor and characteristic 
quantity sets g; obtained correspondingly to the extrinsic factor sets Vj, respectively, each 
consisting of at least one characteristic quantity. The physical model is expressed as a 

15 function f (vj, P) of the extrinsic factor set Vj and a parameter set P consisting of a 
plurality of parameters. An error function S is obtained by multiplying a difference 
between a measured characteristic quantity set g s (s is any one number of i) and a 
calculated function f (v s , P) by a weighting function w s and then summing the squares of 
the products for all extrinsic factor sets v;. Then, a combination of parameters that 

20 minimizes the value of the error function S is extracted. 

Taking an MIS transistor as an exemplary non-linear element, when the gate 
length is longer than about 1 n m, for example, a Frohman-Bentchkowsky model can be 
used as a physical model. Examples of the extrinsic factors include an operating 
temperature x and potentials of a gate electrode (gate voltage V gs ) and a drain electrode 

25 (drain voltage V ds ) with respect to a source electrode. Examples of the characteristic 
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quantities include a current (drain current Ids) flowing between the source electrode and 
the drain electrode and conductance d I ds / d V ds obtained by differentiating the drain 
current Ids by the drain voltage V ds . Examples of the parameters include a threshold 
voltage Vth and a factor j3 as described later. Thus, the number of extrinsic factors 
5 constituting the extrinsic factor set, the number of characteristic quantities constituting 
the characteristic quantity set, and the number of parameters constituting the parameter 
set do not generally agree with one another. 

In the operation of the MIS transistor, for example, the drain current Ids in a 
linear region with a small drain voltage V ds is obtained by approximately j3 (V gs - V t h - 

10 Vds/ 2)V ds . The factor i3 is a value obtained by dividing a product of the capacitance 
C ox of gate insulating film per unit area, the mobility fx of carriers, and the channel 
width W by the channel length L ( j3 = C ox [x W/L). Based on such a model, the 
dependency of the drain current I ds on the gate voltage V gs is obtained from the measured 
values, and the threshold voltage V t h and the factor jS are obtained by calculating 

15 extrapolation and gradients. 

The error function S of the MIS transistor as an example of a non-linear 
element is expressed by the following equation (1): 

Each of the extrinsic factor sets vi, v 2 , V3, v m consists of x 1) extrinsic 
20 factors. When x = 3, for example, the operating temperature x , the gate voltage V gs , 
and the drain voltage V ds can be adopted as extrinsic factors. 

The characteristic quantity set g is a physical quantity of the non-linear element 
with the extrinsic factor set Vj. The drain current I ds and the conductance d I ds / d V ds 
can be adopted as characteristic quantities, for example; and in this case, the number of 



characteristic quantities constituting a characteristic quantity set is two and all of gi, f (v t , 
P) and w; are vector quantities. Since the calculation of the equation (1) is carried out 
for each factor of the vectors, however, the error function S corresponds to the scalar 
quantity. 

5 Further, as characteristic quantities, device simulation results can be used 

instead of measured values. In such a case, generally, a characteristic quantity set can be 
also expressed generally as gi = g (Vj), where g is a function approximating devices by 
device simulation. 

The parameter set P consists of n 2) parameters pi, p2, p3, . .., p n - When n 
10 =2, for example, the threshold voltage and the factor £ can be adopted as 
parameters. 

The weighting functions w; are set correspondingly to the extrinsic factor sets v l5 
respectively. By way of exception, the weighting functions Wj can be set to a constant 1 
for any of the extrinsic factor sets Vj and the error function S can be defined as an absolute 

15 error. Further, when Wj = 1/gj, the error function S is defined as a relative error. 

Parameter extraction is such an operation as above to obtain a parameter set P 
that minimizes the error function S or reduces it to zero within a predetermined error limit. 
The extraction operation proceeds by reducing the value of the error function S while 
updating the parameter set P using a predetermined rule. 

20 To find a combination of parameters that minimizes the value of the error 

function S, a descending method such as Newton's method (in the present specification, 
referred to as "Newton-based solution") and a global search algorithm such as Simulated 
Annealing method (in the present specification, referred to as " combinatorial 
optimization method ") have been conventionally adopted. 

25 The above background-art methods of extracting physical model parameters, 



5 

however, have the following problems. 

The first problem is as follows: in searching the minimum value of the error 
function S by using the Newton-based solution or the combinatorial optimization method 
with a characteristic quantity set g; consisting of a plurality of characteristic quantities, 
5 when there is a variation in size among the values of the errors (referred to as "error 
factors") of the characteristic quantities, with respect to a characteristic quantity whose 
error factor has a large value, the error factor can he effectively reduced but with respect 
to a characteristic quantity whose error factor has a small value, the error factor can not be 
effectively reduced. 

10 In the case, such as the above MIS transistor, where the drain current I<j s and the 

conductance d W d Vd S are adopted as the characteristic quantities constituting the 
characteristic quantity set g„ respective error factors of these characteristic quantities are 
obtained and the search of the minimum value of the error function S is performed by 
multiplying the sum of the squares of these error factors by the weighting function. 

15 Since the error value in the conductance d W d Yds is larger than that in the drain 
current Ids as discussed later in detail, however, the calculation proceeds in a tendency to 
reduce the error factor in the conductance d W d Vd S more than that in the drain current 
Ids- As a result, the accuracy in fitting of the drain current Id S becomes lower than that in 
the case where the drain current Ids is adopted alone as the characteristic quantity. 

20 The drain current Ids in a linear region is obtained by approximately j3 (V gs - 

Vth - Vd S / 2) Vds as discussed earlier. In a saturation region (where V ds > V^: Vd sa t is a 
value of drain voltage V ds at the time when the drain current Id S starts not increasing even 
as the drain voltage V ds increases), however, the drain voltage V ds in the above formula is 
replaced by Va sa t and further the channel length L in the factor /3 is replaced by (L - AL) 

25 (AL represents a spread width of a depletion layer at a drain end) in consideration of 
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channel length modulation effect. Furthermore, AL is a function of (Vd s - V dsat ). 

In this case, though the drain current Ids depends on AL, the conductance 
which is obtained by differentiating the drain current Ids by the drain voltage Vd S depends 
on not only A L but also d A LI d V ds . Though it is possible to express the behavior of 
5 the AL relative to the variation in (Va s - Vd sa t) with high accuracy from depletion 
approximation, it is likely to cause a problem in accuracy of expression for d A L/ d Vd S 
especially near Vd s = Vdst- Therefore, the accuracy of the conductance which includes 
d A L/ 3 V ds as a parameter is lower than that of the drain current Ids which is not 
affected by that term. 

10 Specifically, since the model formula in the MIS transistor is set by connecting 

the formulae of the drain currents in the operating regions (weak inversion region, strong 
inversion/linear region, strong inversion/saturation region), it is likely to cause 
deterioration in accuracy of the value of conductance which is obtained by differentiating 
the drain current by the drain voltage, especially near the connections of the regions. 

15 Therefore, the value of error factor in the conductance d W 3 Vd S is likely to becomes 
larger than that in the drain current Id S . As a result, when fittings of these error factors 
are performed simultaneously, there is a tendency to reduce the error factor in the 
conductance 3 Ids/ 3 Vd S more than that in the drain current Id S . 

On the other hand, the second problem in the background-art method of 

20 extracting physical model parameters is as follows: when the combinatorial optimization 
method is used in searching the minimum value of the error function S, its computational 
expenses for reaching a true minimum value disadvantageous^ become enormous, 
though it is possible to avoid entrapment in local solutions. Further, use of the 
Newton-based solution causes entrapment in local solutions. 



25 



7 

SUMMARY OF THE INVENTION 

The present invention is directed to a method of extracting physical model 
parameters. According to a first aspect of the present invention, the method of 
extracting physical model parameters comprises the steps of: (a) adopting a physical 
5 model for physical properties that provide characteristic quantity sets gi (i = 1, 2, m) 
each consisting of the first to z-th (z ^ 2) characteristic quantities gj y (y = 1, 2, z), 
corresponding to a plurality of extrinsic factor sets Vj, respectively, each of which consists 
of at least one extrinsic factor, the physical model expressing a calculated value of each 
characteristic quantity set g s (s is any one value of i) as a function f y (v s , P) of 

10 corresponding one extrinsic factor set v s and a parameter set P consisting of a plurality of 
parameters; and (b) obtaining an error function S by summing values each of which are 
obtained by dividing the square of the difference between each of the characteristic 
quantities gj y and corresponding function f y (v;, P) by variance Oj y 2 of observed values of 
the characteristic quantities gi y obtained by observing the physical properties of a plurality 

15 of samples, for the plural number of extrinsic factor sets and extracting the parameter set 
P that gives the minimum value of the error function S. 

According to a second aspect of the present invention, in the method of 
extracting physical model parameters of the first aspect, the step (b) is the step of: (b-a) 
extracting the parameter set P that gives the minimum value of the error function S by 

20 obtaining the error function S repeatedly with the parameter set P updated, and the step 
(b) includes the step of: (b-1) verifying if the function f y (v s , P) obtained by calculation 
using updated parameter set P reproduces the observed values of each characteristic 
quantity set g s obtained from the plurality of samples by utilizing conformity of the value 
of the error function S to x 2 distribution and extracting the parameter set P at that time as 

25 one that gives the minimum value of the error function S when the function f y (v s , P) 
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reproduces the observed values. 

Preferably, the step (b) is the step of: obtaining the error function S repeatedly 
with the parameter set P updated while maintaining a probability Q of increasing the 
value of the error function S positive. 
5 Preferably, the probability Q is obtained on the basis of a predetermined 

amount which fluctuates monotonously as the parameter set P is updated, and the 
predetermined amount contributes to a tendency to decrease the probability Q as the 
parameter set P is updated. 

Preferably, the step (b) is the step of: updating the parameter set P only in such 
10 a direction that the error function S decreases monotonously. 

Preferably, Gauss-Newton method is adopted in the step (b). 

Preferably, Levenberg-Marquardt method is adopted in the step (b). 

According to a third aspect of the present invention, in the method of extracting 
physical model parameters of the second aspect, the step (b) further includes the steps of: 
15 (b-2) judging if the value of the error function S obtained with the parameter set P 
updated is considered as the minimum value when the function f y (v s , P) does not 
reproduce the observed values of the characteristic quantity set g s obtained from the 
plurality of samples in the step (b-1) and extracting the parameter set P at that time as one 
that gives the minimum value of the error function S when the value of the error function 
20 S is considered as the minimum value; and (b-3) judging if the number of updates 
exceeds a predetermined number of times when the value of the error function S is not 
considered as the minimum value in the step (b-2) and extracting the parameter set P at 
that time as one that gives the minimum value of the error function S when the number of 
updates exceeds the predetermined number of times or obtaining the value of the error 
25 function S with the parameter set P updated to go back to the step (b-1) when the number 
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of updates does not exceed said predetermined number of times. 

According to a fourth aspect of the present invention, in the method of 
extracting physical model parameters of the third aspect, the method used in the step (b-a) 
of extracting the parameter set P that gives the minimum value of the error function S is 
5 changed to execute the steps (b-1) to (b-3) again, instead of extracting the parameter set P 
at that time as one that gives the minimum value of the error function S, when it is judged 
that the number of updates exceeds the predetermined number of times in the step (b-3). 

The present invention is also directed to a computer-readable storage medium. 
According to a fifth aspect of the present invention, the computer-readable storage 

10 medium stores a program that causes a computer to execute the method of extracting 
physical model parameters as defined in any one of the first to fourth aspects 
independently or in combination with a program previously stored in the computer. 

The present invention is still directed to a method of manufacturing a non-linear 
element. According to a sixth aspect of the present invention, the method of 

15 manufacturing a non-linear element manufactures a non-linear element by performing: a 
characteristic simulation which adopts device modeling using the method of extracting 
physical model parameters as defined in any one of the first to fourth aspects; and a 
physical process on the basis of the characteristic simulation. 

In the method of extracting physical model parameters of the first aspect, since 

20 the error function S is obtained by summing values each of which are obtained by 
dividing the square of the difference between each of the characteristic quantities gi y and 
corresponding function f y (v;, P) by variance Oj y 2 of observed values from a plurality of 
samples for the plural number of extrinsic factor sets in the step (b), it is possible to 
correct the effects of variation and bias of the observed values included in the error factor 

25 with division by the variance Oj y 2 for each of the characteristic quantities even when there 
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is a variation in size among the values of the error factors in the characteristic quantities 
in the error function S. Therefore, the error function S becomes a function normalized 
on each characteristic quantity, and that prevents a strong effect on the error function S by 
biased error factor of a specified characteristic quantity and provides the method of 
extracting physical model parameters which allows sufficient reduction in error value for 
each characteristic quantity. 

In the method of extracting physical model parameters of the second aspect, 
since it is verified whether the function f y (v s , P) reproduces the observed values of each 
characteristic quantity set g s by utilizing conformity of the value of the error function S to 
X 2 distribution in the step (b-1) and the parameter set P at that time is extracted as one that 
gives the minimum value of the error function S when it is verified that the function f y (v s , 
P) reproduces the observed values, it is not necessary to repeat the calculation until the 
value of error function S converges to be considered as minimum and therefore an 
efficient and quick extraction of the parameter set P can be achieved. 

In the method of extracting physical model parameters of the third aspect, since 
the step (b) includes the step (b-2) of judging if the value of the error function S obtained 
with the parameter set P updated in the step (b-1) is considered as the minimum value and 
the step (b-3) of judging if the number of updates exceeds a predetermined number of 
times, a subsidiary extraction of the parameter set P can be achieved even when it is 
verified that the function f y (v s , P) does not reproduce the observed values. 

In the method of extracting physical model parameters of the fourth aspect, 
when it is judged that the number of updates exceeds the predetermined number of times 
in the step (b-3), the method used in the step (b-a) is changed to execute the steps (b-1) to 
(b-3) again, instead of extracting the parameter set P at that time as one that gives the 
minimum value of the error function S. Therefore, it becomes possible to adopt one of 
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the combinatorial optimization method and the Newton-based solution first in the step 
(b-a) and adopt the other one in the step (b-a) to retry the extraction of the parameter set P 
when the first method fails the extraction of the parameter set P. 

In the computer-readable storage medium of the fifth aspect, it is possible to 
cause a computer to execute the method of extracting physical model parameters as 
defined in any one of the first to fourth aspects. 

In the method of manufacturing a non-linear element of the sixth aspect, since 
the physical process are carried out on the basis of the characteristic simulation with high 
accuracy at low computational cost, the non-linear element can be manufactured, being 
close to a design specification, at low cost. 

An object of the present invention is to provide a method of extracting physical 
model parameters, which allows sufficient reduction in value of an error factor for each 
characteristic quantity even if there is a variation in size among the values of the error 
factors of the characteristic quantities in the error function, and further provide a 
technique to perform an efficient and quick parameter search, i.e., parameter extraction 
for obtaining a true solution. 

Another object of the present invention is to provide a storage medium for 
storing a computer-driven program to extract parameters. 

Still another object of the present invention is to provide a method of 
manufacturing a non-linear element, which includes physical simulations using the 
parameter extraction technique. 

These and other objects, features, aspects and advantages of the present 
invention will become more apparent from the following detailed description of the 
present invention when taken in conjunction with the accompanying drawings. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

Fig. 1 is a flowchart showing an overview of manufacturing process of an LSI 
device in accordance with the present invention; 

Fig. 2 is a flowchart showing an operation of a parameter extractor in 
accordance with the present invention; 

Fig. 3 is a flowchart showing a method of extracting physical model parameters 
in accordance with a first preferred embodiment of the present invention; 

Fig. 4 is a flowchart showing a Newton-based solution; 

Fig. 5 is a flowchart showing a method of extracting physical model parameters 
in accordance with a second preferred embodiment of the present invention; and 

Figs. 6 and 7 are flowcharts showing a method of extracting physical model 
parameters in accordance with a third preferred embodiment of the present invention. 

DESCRIPTION OF THE PREFERRED EMBODIMENTS 
< The First Preferred Embodiment > 
A. Basic Idea: 

The first preferred embodiment shows a method of extracting physical model 
parameters, which allows sufficient reduction in value of an error factor for each 
characteristic quantity even if there is a variation in size among the values of the error 
factors of the characteristic quantities in the error function. 

Prior to giving a detailed description of the first preferred embodiment, we will 
describe a basic idea of the present invention. Also in the first preferred embodiment, 
like in the background-art technique, the physical model is set as the physical properties 
of a non-linear element, where a non-linear relation is established between extrinsic factor 
sets Vj (i = 1, 2, 3, m) each consisting of at least one extrinsic factor and 
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corresponding characteristic quantity sets gj each consisting of at least one characteristic 
quantity. The physical model is expressed as a function f (v is P) of the extrinsic factor 
set V; and a parameter set P consisting of a plurality of parameters. Further, an error 
function S is obtained on the basis of this physical model. 

In the first preferred embodiment, however, the sum of values obtained by 
dividing the squares of the error factors which are obtained by observation or calculation 
in device simulation of the properties of u (u ^ 2) prototypes by variances of observed 
values of corresponding characteristic quantities for the plural number of extrinsic factor 
sets (e.g., total number m) is adopted as the error function S of the physical model. 

Then, a combination of parameters that minimizes the value of the error 
function S is extracted. 

The present invention is not to be restricted to the application to the 
semiconductor field but is also applicable to other fields such as electricity, machinery, 
and chemistry if it is the technique for adopting a physical model expressed as a function 
containing parameters as above described. The following discussion will be made, 
taking a field of method of manufacturing semiconductor devices as an example. 

B. Application to Method of Manufacturing Semiconductor Device: 
bl) Overview of Method of Manufacturing Semiconductor Device. 
Fig. 1 is a flowchart showing an overview of the manufacturing process of an 
LSI device to which the present invention is applicable. The manufacturing process is 
broadly divided into a design process group 90 and a semiconductor processing step 905 
which is a physical process. The design process group 90 is broadly divided into a 
functional design process 901, a logic design process 902, a circuit design process 903, 
and a layout design process 904. In the semiconductor processing step 905, a 
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semiconductor processing is performed on the basis of information from the design 
process group 90, and an LSI device 300 is thereby obtained. 

The execution of the circuit design process 903 is carried out by a circuit 
simulator 1 and a parameter extractor 3. In addition, a timing simulator 201 may be 
adopted. The circuit simulator 1 is a main device for performing circuit simulations, and 
the parameter extractor 3, which supplies parameters to the circuit simulator 1, receives 
measured values of the manufactured LSI device 300 and simulation results from a device 
simulator 202. The circuit design process 903 may also adopt a parameter database 2 for 
storing parameters determined by the parameter extractor 3. In the example of Fig. 1, 
the circuit design process 903 adopts the circuit simulator 1, the parameter database 2, the 
parameter extractor 3, and the timing simulator 201, all of which are enclosed by a box 
indicating the circuit design process 903. 

Fig. 1 is merely one example in schematic form; so the process steps 901 to 904 
in the design process group 90 need not be independent of one another and the parameter 
extractor 3 needs not be individual hardware. For example, the whole design process 
group 90 may be implemented by a single computer operating on the basis of a 
predetermined program. It is, of course, possible to operate the parameter extractor 3 by 
exclusive software or by a patch program for existing software to operate the whole 
design process group 90. Such software and program can be stored in a 
computer-readable storage medium. 

Fig. 2 is a flowchart showing an operation of the parameter extractor 3 in 
accordance with the present invention. First, measured values obtained by measuring 
equipment such as a parameter analyzer and an LCR meter are inputted to an initial-value 
determining unit 11. It is preferable that the measuring equipment should operate under 
automatic control. Alternatively, device simulation results obtained by the device 
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simulator 202 may be inputted. 

To determine parameters in device modeling, repeated calculations are usually 
performed. For efficient repeated calculations, the initial-value determining unit 11 
determines initial values of such parameters. A proper selection of these initial values is 
deeply involved in the precision of parameters to be finally obtained. For this reason, 
the initial-value determining unit 1 1 adopts a simple model which restricts the operating 
range and shape of devices to explicitly determine the initial values of some parameters 
without necessitating the repeated calculations. 

In the operation of an MIS transistor, for example, the drain current I ds in a 
linear region with a small drain voltage V ds is obtained by approximately j3 (V gs - V t h - 
Vds/ 2)Vd s . Based on such a model, the dependency of the drain current I ds on the gate 
voltage V gs is obtained from the measured values, and the threshold voltage V th and the 
factor 0 are obtained by calculating extrapolation and gradients. 

The initial values of some parameters obtained in this way are given to a 
parameter optimization unit 12 together with the measured values or device simulation 
results, where parameters are determined. The details of an operation of the parameter 
optimization unit 12 will be discussed later. 

Device characteristics calculated by using the parameters obtained in the 
parameter optimization unit 12 are displayed in a precision verification unit (display 
device) 13 together with those obtained from the measured values, such as the 
dependency of the drain current I ds on the drain voltage V ds . This allows visual 
recognition of the precision of determined parameters. In Fig. 2, the dependency of the 
drain current I ds on the drain voltage V ds is plotted at different gate voltages V u V 2 and 
V 3 . 

When parameters with satisfactory precision are acknowledged, the parameters 
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are stored in the parameter database 2 for reuse. 

Fig. 2 is merely one example in schematic form; so the units 1 1 to 13 need not 
be independent of one another. For example, the whole parameter extractor 3 may be 
implemented by a single computer operating on the basis of a predetermined program. It 
is, of course, possible to operate the parameter optimization unit 12 by exclusive software 
or by a patch program for existing software to operate the whole design process group 90. 
Such software and program can be stored in a computer-readable storage medium. That 
allows a computer to execute the method of extracting physical model parameters. 

Further, when a characteristic simulation is performed by using the method of 
extracting physical model parameters of the first preferred embodiment and the physical 
process is executed on the basis of the characteristic simulation to manufacture a 
non-linear element such as an MIS transistor, since the physical process is executed on 
the basis of the characteristic simulation with a high degree of accuracy and low 
calculation cost, it is possible to manufacture a non-linear element close to design 
specifications at low cost. 

b2) Overview of Parameter Extraction. 

Fig. 3 is a flowchart showing an operation of the parameter optimization unit 12 
in an exploded form. This flowchart shows a step of extracting the parameter set P that 
gives the minimum value of the error function S described in the section A. 

In the step SOI, a parameter search is performed by a minimum-value solver 
such as the combinatorial optimization method or the Newton-based solution. Herein, 
the minimum-value solver refers to a solution for extracting the parameter set P that gives 
the minimum value of the error function S by obtaining the error function S repeatedly 
with the parameter sets P updated. As an example of the minimum-value solver, the 
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combinatorial optimization method will be discussed in the section b3) and the 
Newton-based solution will be discussed in the section b4) in detail. 

Then, in the step S02, a judgment is made on whether the value of the error 
function S is considered as the minimum value or not. The judgment may be performed 
5 according to a convergence condition for each minimum-value solver. When it is 
judged that the value of the error function S is sufficiently small, considering that the 
error function S converges, the parameter set P at that time is extracted as one that gives 
the minimum value of the error function S (step S04). 

On the other hand, when the value of the error function S is not considered as 
10 the minimum value, a judgment is made in the step S03 on whether the number of 
updates of the parameter set P, #Iter (its initial value is 0), exceeds a predetermined 
number of times Imax or not. When #Iter exceeds Imax, considering that the error 
function S converges, the parameter set P at that time is extracted as one that gives the 
minimum value of the error function S. On the other hand, when #Iter does not exceed 
15 Imax, the process returns to the step SOI and the error function S is obtained again with 
the parameter set P updated by using the minimum-value solver. 

Further, "#Iter++" in the flowchart of Fig. 3 represents increment of the number 
of updates #Iter. 

In the first preferred embodiment, as discussed in the section A, the error 
20 function S of the non-linear element such as the MIS transistor is expressed by the 
following equation (2): 




(/ = 1, 2, —,m) 
(y = l,2,-,z) 



where the extrinsic factor set v ; , the characteristic quantity g iy , the parameter set P and the 
function f y (v ; , P) are the same as those in the background art, and the variance a iy 2 is a 
variance of the observed values of the y-th characteristic quantity in the characteristic 
quantity set gi. Herein, the subscript y is exhibited in order to explicitly indicate that 
there are y (y = 1, 2, z) kinds of characteristic quantities, and the equation expresses 
that summation for the number y is made. 

The s-th (s is any one value of i) characteristic quantity set g s is expressed by 
the following equation (3): 

X s m (y s ) 



\g S2 

(3) 

(q = l,2,-,u) 

where g syq (v s ) represents an observed value of the y-th characteristic quantity of 
the characteristic quantity set g s obtained when the extrinsic factor set v s is given to the 
q-th sample and u denotes the number of samples. Taking a MIS transistor which is 
singly formed on a chip as an example, the observed value of the first characteristic 
quantity, the drain current I ds , in the q-th sample out of a plurality of chips obtained by 
actual device manufacture or device simulation under the condition of the extrinsic factor 
set v s , for example, corresponds to g s i q (v s ). Further, as can be understood from the 
equation (3), the characteristic quantity g sy is an average value of the observed values 
gs yq (v s ) of each physical property of the samples. 

Furthermore, the variance o s 2 of the characteristic quantity set g s and the 
variance a sy 2 of each characteristic quantity are expressed by the following equations 
(4): 
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Z {gsyg (V s )-g sy } 2 




Herein, the variance a sy 2 is an unbiased estimation value. Therefore, the 
equation (4) is a value obtained by multiplying the sample variance by u/(u-l). 
Specifically, though u is used in the denominator of the equation (4), instead of (u-1), in 
the case of sample variance, the value obtained by multiplying the sample variance by 
u/(u-l) is adopted in the equation (4) in order to express the variance of the population of 
the samples as an unbiased estimation value. 

The reason why the variance a \ y 2 of the observed values of the characteristic 
quantities g iy in the characteristic quantity set g, is introduced in the error function S is as 
follows. 

As indicated in the equation (2), by dividing the square of error factor of each 
of the characteristic quantities g iy by the variance o iy 2 of the observed values of a 
plurality of samples, it is possible to correct the effects of variation and bias of the 
observed values included in the error factor with division by the variance a iy 2 for each of 
the characteristic quantities gi y even when there is a variation in size among the values of 
the error factors in the characteristic quantities g iy . In other words, since the error 
function S is a function normalized on each of the characteristic quantities g iy , it is 
possible to prevent a strong effect on the error function S by biased error factor in specific 
one of the characteristic quantities g, y . 

Therefore, it is possible to achieve the method of extracting physical model 
parameters, which does not cause the first problem of the background art even when there 
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is a variation in size among the values of the error factors in the characteristic quantities 
and allows sufficient reduction in value of the error factor for each characteristic quantity. 

Further, the process described in the steps SOI to S04 of Fig. 3 may be 
performed on a computer by running a program for integrally processing the steps. 
Alternatively, the process described in the steps SOI to S04 may be integrally performed 
by using hardware. The process may be carried out on a computer by running a program 
for executing the steps SOI to S04 separately. 

Alternatively, hardware executing the steps SOI to S04 separately may be used 
as the parameter optimization unit 12. 

b3) Combinatorial Optimization Method. 

Examples of the combinatorial optimization method include well-known 
techniques called simulated annealing and simulated diffusion (hereinafter referred to as a 
"SA method" and "SD method", respectively). The application of the SA method to 
semiconductor devices is discussed, for example, by Man-Kuan. Vai, et al in an article 
entitled "Modeling of Microwave Semiconductor Devices Using Simulated Annealing 
Optimization", IEEE Trans. Electron Devices, vol. ED-36, No. 4, pp. 761-762, Apr. 1989 
(hereinafter referred to as a "document 1"). The application of the SD method to 
semiconductor devices is discussed, for example, by T. Sakurai, et al. in an article entitled 
"Fast Simulated Diffusion: An Optimization Algorithm for Multiminimum Problems and 
Its Application to MOSFET Model Parameter Extraction", IEEE Trans. Computer- Aided 
Design, vol. CAD-11, No. 2, pp. 228-233, Feb. 1992 (hereinafter referred to as a 
"document 2"). 

In the SA method and the SD method, the error function S is obtained from 
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parameters incremented by a predetermined amount, and when the error function S is 
increased, the incremented parameters are adopted as updated parameters with a 
probability Q. In the document 1, for example, an amount for updating parameters 
(hereinafter referred to as the "amount of update") is expressed as A V = RV 0 , where V 0 
is the basic amount of update and R (0 ^R ^ 1) is a random number. The parameters 
updated by A V are used for the next calculation unless they are beyond the specified 
limits. In the document 2, the amount of update proportional to the gradient of the error 
function S is added to parameters. When the error function S obtained from resultant 
parameters is increased, these parameters are adopted as updated parameters with a 
probability Q. More specifically, when the obtained parameters increase the value of the 
error function S, a random number R is generated, and these parameters are adopted as 
updated parameters on condition that the value of random number is lower than the 
probability Q. When the value of the error function S decreases, on the other hand, 
parameter update is performed. 

Both the SA method and the SD method adopt a so-called Metropolis algorithm. 
The Metropolis algorithm is introduced for example by R. A. Rutenbar in an article 
entitled "Simulated Annealing Algorithms: An Overview", IEEE Circuits and Devices 
Magazine, pp. 19-25, Jan. 1989 (hereinafter referred to as a "document 3"). Referring to 
this algorithm along with the SA method and the SD method, a probability Q that the 
error function S increases by the amount of increase 5 S (> 0) is expressed by exp(- 8 
S/T), where 0 < 1 holds. 

The divisor T is a predetermined amount that decreases as the calculation is 
repeated with an update of the parameters. In the document 1, for example, it is called a 
pseudo-temperature, the initial value of which is set to 500 or more and updated to 90 
percent of the value as the repeated calculation is performed. 
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A convergence condition adopted in the SA method and the SD method is, for 
example, whether the magnitude of parameters is within the specified limits, whether the 
number of times random number generation occurs reaches a predetermined upper limit, 
or whether the pseudo-temperature T is cooled enough. 

Additionally, as a convergence condition of the combinatorial optimization 
method, whether the error function S decreases a predetermined number of times in 
succession may be adopted. 

Specifically, the combinatorial optimization method should be terminated when 
the error function S has decreased a predetermined number of times in succession, in 
which case it is judged that the error function S has decreased monotonously with each 
calculation carried out with the update of parameters. That is, the convergence 
condition can be expressed by the following equation (5): 

Sk > Sk+i > . . . Sk+t (5) 
where S k is the value of the error function S obtained with k repeated calculations. 

The number of times t that the value of the error function S has decreased in 
succession is, for example, set to 6. This convergence condition corresponds to the step 
S02 of Fig. 3. When the error function S has decreased during t repeated calculations 
carried out with the update of the parameter set P, the process goes to the step S04; 
otherwise, the process goes to the step S03. 

Besides the Metropolis algorithm, the SD method further adopts the term of the 
Brownian movement as part of the amount of update of parameters. In the document 2, 
for example, the amount of update dx is expressed by the following equation (6): 
dx = -Vf(x)dt + V2Tdw (6) 

In the equation (6), the parameter x and the function f are equivalent to the 
parameter set P and the error function S, respectively. In the right side, the first term is 
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the drift term and the second term is equivalent to the Brownian movement, 
is the Gaussian random noise. 



Herein, dw 



b4) Newton-based Solution. 



5 



The following equation (7) should be satisfied to minimize the error function S: 



— = 0 (j = l,2,...,n) 



This is n-order non-linear simultaneous equations, where the parameters pi, p2, 
p n constituting the parameter set P are unknown variables. As a technique for 
solving the simultaneous equations, Newton's method is the most common numerical 
10 solution, in which the parameter set P is updated after each calculation so that the error 
function S obtained from the updated parameter set P decreases monotonously. In this 
solution, the amount of update A P with respect to the parameter set P is calculated by 
using the following equation (8): 



AP= [j T WJ - H (k) W(G - F)Y J T W{G - F) 




where 




_ d 2 /(v„P) 
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Further, k = 1, 2, n and J T is the transposed matrix of J. 



The terms inside 



the square brackets indicate the Hessian of the function f (v, P). 



To improve the efficiency of calculation of the amount of update AP , various 
approximation methods are applied. In the Gauss-Newton method, for example, the 
terms corresponding to second and higher differentiations of the function f are discarded 
as follows: 

where: A = W 1/2 J, B = W I/2 (G-F), C = A T B 

According to the equation (9), the parameter set P is sequentially obtained by 
QR decomposition, e.g. by Householder method. 

For further improvement in convergence, the Levenberg-Marquardt method has 
been proposed, wherein the diagonal terms are added as approximation of the Hessian to 
10 emphasize an element in such a direction that the error function S decreases. This 
technique is discussed for example by K. Dogains and D. L. Scharfetter in an article 
entitled "General Optimization and Extraction of IC Device Model Parameters", IEEE 
Trans. Electron Devices, vol. ED-30, No. 9, pp. 1219-1228, September 1983. More 
specifically, the amount of update AP is calculated on the basis of the following 
15 equation (10): 

ap = [a t a + ^d]- 1 c (iq) 

where D = I + diag(A T A) 
where I is an unit matrix and diag() represents a matrix that adopts diagonal elements of a 
matrix in the parentheses as diagonal elements and zero as the other elements. The 
factor X is set to a large value (e.g., approximately 0.1) at the beginning of the repeated 
20 calculations and is set to zero near the solution. 

Fig. 4 is a flowchart of the Newton-based solution using the equation (9). In 
the step 124, the parameter set P (r) obtained in the r-th calculation is increased by the 
amount of update {[A^J^C} to calculate a parameter set P (r+1) to be obtained in the 
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(r+l)th calculation. In the step 125, the error function S is updated using the parameter 
set ~P^ T+1 \ In the step 126, it is determined whether the error function S is small enough, 
i.e., whether the error function S is zero within the specified limits. When the error 
function S becomes zero within the specified limits, the optimization process is 
5 terminated; otherwise, the process returns to the step 124 after the number of repetitions r 
is incremented by one in the step 127. 

If it is judged in the step 126 that the updated parameter set P (r+1) and the 
parameter set P (r) before update are within the specified limits, the optimization process is 
terminated since further repetition of calculations would not improve the precision of the 
10 parameter set P. This is a convergence condition for the Newton-based solution. If 
they are not within the specified limits, the process returns to the step 124 through the 
step 127. 

The steps 124 and 125 in Fig. 4 correspond to the steps SOI in Fig. 3 and the 
step 126 in Fig. 4 corresponds to the step S02 in Fig. 3. 

15 

C. Variation: 

While the document 3 discusses the case of adopting a Boltzmann distribution 
to show the probability in the Metropolis algorithm, a Lorentz distribution can be adopted 
as discussed in the document 2. Alternatively, it is also possible to adopt a function 
20 which discards terms of higher order than a term of a predetermined order in a series 
obtained by expanding the Boltzmann distribution. 
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< The Second Preferred Embodiment > 

The second predetermined embodiment is a variation of the method of 
extracting physical model parameters in accordance with the first preferred embodiment, 



26 

and is a method of extracting physical model parameters which allows an efficient and 
quick parameter extraction as well as sufficient reduction in value of an error function for 
each characteristic quantity. 

When fitting of the function f y (vj, P) set as a physical model and an observed 
5 value is performed by using the error function S shown in the equation (2), there is a 
problem to what degree the value of the error function S should be reduced for the 
judgment that a high-quality parameter set P can be obtained. 

Specifically, as the second problem discussed in the description of the 
background art, there is a problem that computational expenses for reaching a true 

10 minimum value become enormous when the combinatorial optimization method is used 
as the minimum-value solver. Further, there is another problem that use of the 
Newton-based solution causes entrapment in local solutions. 

Then, when it is judged by some standard that minimization of the error is 
sufficiently made, it is efficient that the calculation is immediately terminated at that stage 

15 and the parameter set P is extracted. 

When attention is paid to the error function S shown in the equation (2), 
assuming that the function f y (vi, P) represents the average value of the characteristic 
quantities g; y with respect to a population of u samples to be observed, it is considered 
that the equation (2) conforms to ^ distribution. The reason is as follows: assuming that 

20 the error factor of the observed value of the characteristic quantity g sy is incidental at each 
bias point (e.g., v s ) of the extrinsic factor sets Vj, it is considered that the error distribution 
in the population is a normal distribution, and the characteristic quantity g sy in the 
equation (2) represents the observed value of the physical property of each sample and the 
variance a sy 2 represents an unbiased estimation value of the variance of the population. 

25 Therefore, utilizing the conformity of the error function S shown in the 
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equation (2) to % 2 distribution, it is verified if the function f y (v s , P) obtained by 
calculation using the updated parameter set P reproduces the observed values of the 
characteristic quantity set g s obtained from a plurality of samples, and when it is verified 
that the function reproduces it, the parameter set P at that time is extracted as one that 
5 gives the minimum value of the error function S. 

For this verification, a judgment is made on whether the hypothesis that the 
function f y (v;, P) represents the average value of the population to be observed is rejected 
or not. Specifically, the judgment on whether the hypothesis is rejected or not is made 
depending on whether the value of probability obtained by integrating % 2 distribution 

10 function from 0 to the value of the error function S calculated on the basis of the value of 
the updated parameter set P falls within a predetermined standard for rejection. When 
the hypothesis is rejected, update of the parameter set P is performed, instead of adopting 
the parameter set P at that time. On the other hand, when the hypothesis is not rejected, 
considering that the function f y (v s , P) reproduces the observed values of the characteristic 

15 quantity set g s , it is judged that the parameter set P at that time has high quality. 

More specifically, the verification is performed according to the following 
procedure. First, the % 2 distribution function is expressed by the equation (11): 

(-) (f_1) 

F X (S) = -Z exp(-|) (11) 

2T(|) 2 

where F X (S), x and T represent % 2 distribution function, degree of freedom and 
20 gamma function, respectively. Further, the degree of freedom x is calculated by using 
the number m of bias points of the extrinsic factor sets Vi and the number n of parameter 
sets P, being expressed as x = m - n - 1 . 

Then, by integrating the equation (11) from S = 0 to S = S 0 (S 0 is a value of 



error function S obtained by calculation on the basis of the values of the updated 
parameter set P) as shown in the left side of the following equation (12), a probability 
value of % 2 distribution function is calculated. 

^F x (S)ds<(l-a) (12) 

5 where (1 - a ) of the right side represents the rejection standard, and as a , values such 
as 0.01, 0.05 and 0.10 are adopted. 

As shown in the equation (12), when the calculated value in the left side 
becomes smaller than the value in the right side, the hypothesis is not rejected and the 
parameter set P at that time is adopted as one that gives the minimum value of the error 
10 function S. On the other hand, when the calculated value in the left side becomes larger 
than the value in the right side, the hypothesis is rejected, not adopting the parameter set 
P at that time, and it is verified that the function f y (v s , P) does not reproduce the observed 
values of the characteristic quantity set g s . 

Further, when it is verified that the function f y (v s , P) obtained by calculation 
15 using the updated parameter set P does not reproduce the observed values of the 
characteristic quantity set g s , a judgment is made on whether the parameter set P at that 
time gives the minimum value of the error function S according to the convergence 
condition for the minimum- value solver like in the first preferred embodiment. 

Then, when it is judged to give the minimum value, the parameter set P at that 
20 time is adopted as one that gives the minimum value of the error function S and when it is 
not judged to give the minimum value, a judgment is made on whether the parameter set 
P should be updated according to the number of updates like in the first preferred 
embodiment. 

The flowchart of Fig. 5 shows a series of steps as discussed above. 



29 

Specifically, the step Sll is a step for performing the minimum-value solver such as the 
combinatorial optimization method or the Newton-based solution, and the step S 12 is a 
step for calculating the probability value of the x 2 distribution function and verifying if 
the function f y (v s , P) obtained by calculation using the updated parameter set P 
5 reproduces the observed values of the characteristic quantity set g s . The step S13 is a 
step for judging if the convergence condition for the minimum-value solver is satisfied 
and the step S 14 is a step for judging if the parameter set P should be updated according 
to the number of updates. 

When it is verified that the function reproduces the observed values in the step 

10 SI 2 or it is judged that the convergence condition is satisfied in the step S13, considering 
that the error function S converges, the parameter set P at that time is extracted as one 
that gives the minimum value of the error function S (step SI 5). When it is judged that 
the number of updates #Iter does not exceed the predetermined number of times Imax in 
the step SI 4, the parameter set P is updated and the process returns to the step Sll, and 

15 when it is judged that #Iter exceeds Imax, the parameter set P at that time is extracted as 
one that gives the minimum value of the error function S (step SI 5). 

In the method of extracting physical model parameters of the second preferred 
embodiment, since it is verified whether the function f y (v s , P) reproduces the observed 
values by utilizing the conformity of the value of the error function S to % 2 distribution 

20 function and the parameter set P at that time is extracted as one that gives the minimum 
value of the error function S when it is verified that the function reproduces the observed 
values, it is not necessary to repeat the calculation until the value of error function S 
converges to be considered as minimum and therefore an efficient and quick extraction of 
the parameter set P can be achieved. 

25 Further, since the process comprises the steps SI 3 and SI 4, even when it is 
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verified that the function f y (v s , P) does not reproduce the observed values, a subsidiary 
extraction of the parameter set P can be achieved. 



< The Third Preferred Embodiment > 
5 The third preferred embodiment is a variation of the method of extracting 

physical model parameters in accordance with the second preferred embodiment. 
Specifically, in the third preferred embodiment, when the number of updates #Iter 
exceeds the predetermined number of times Imax in the step S14 of the flowchart shown 
in Fig. 5, the method of minimum- value solver is changed, instead of extracting the 

10 parameter set P at that time as one that gives the minimum value of the error function S. 
Then, using a new minimum-value solver, verification on whether the function f y (v s , P) 
reproduces the observed values, judgment on whether the convergence condition for the 
minimum-value solver is satisfied and judgment on whether the parameter set P should be 
updated according to the number of updates are performed again. 

15 Fig. 6 is a flowchart showing a case, as an example of the third preferred 

embodiment, where the combinatorial optimization method is adopted as the 
minimum-value solver first and when the number of updates #Iter exceeds the 
predetermined number of times Imax, the Newton-based solution is adopted as the 
minimum- value solver. Specifically, the step S21 is a step for performing the 

20 combinatorial optimization method and the step S22 is a step for calculating the 
probability value of x 2 distribution function and verifying if the function f y (v s , P) 
obtained by calculation using the updated parameter set P reproduces the observed values 
of the characteristic quantity set g s . The step S23 is a step for judging if the convergence 
condition for the minimum- value solver is satisfied and the step S24 is a step forjudging 

25 if the parameter set P should be updated according to the number of updates. When it is 
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verified that the function reproduces the observed values in the step S22 or it is judged 
that the convergence condition is satisfied in the step S23, considering that the error 
function S converges, the parameter set P at that time is extracted as one that gives the 
minimum value of the error function S (step S29). When it is judged that the number of 
5 updates #Iter does not exceed the predetermined number of times Imax in the step S24, 
the parameter set P is updated and the process returns to the step S21 . 

On the other hand, when it is judged that the number of updates #Iter exceeds 
the predetermined number of times Imax in the step S24, the number of updates #Iter is 
reset to zero and the minimum-value solver is changed to the Newton-based solution 

10 (step S25). Then, the steps S26 to S28 like the steps S22 to S24 are performed. 

When it is verified that the function reproduces the observed values in the step 
S26 or it is judged that the convergence condition is satisfied in the step S27, considering 
that the error function S converges, the parameter set P at that time is extracted as one 
that gives the minimum value of the error function S (step S29). When it is judged that 

15 the number of updates #Iter does not exceed the predetermined number of times Imax in 
the step S28, the parameter set P is updated and the process returns to the step S25 and 
when it is judged that #Iter exceeds Imax, the parameter set P at that time is extracted as 
one that gives the minimum value of the error function S (step S29). 

Fig. 7 is a flowchart showing a case, as another example of the third preferred 

20 embodiment, where the Newton-based solution is adopted as the minimum-value solver 
first and when the number of updates #Iter exceeds the predetermined number of times 
Imax, the combinatorial optimization method is adopted as the minimum-value solver. 
The steps S31 to S39 in Fig. 7 are the same as the steps in Fig. 6 except that the step S21 
for performing the combinatorial optimization method and the step S25 for performing 

25 the Newton-based solution are exchanged. 
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In the case of Fig. 6, since the combinatorial optimization method is performed 
first, the error function S can be globally approximated to the minimum value to some 
degree and after that, the local minimum value can be obtained by performing the 
Newton-based solution. 

5 On the other hand, in the case of Fig. 7, since the Newton-based solution is 

performed first, the computational expenses are small and high-quality parameter set P 
can be sometimes obtained in early stage. Then, even if no high-quality parameter set P 
is obtained by the Newton-based solution, it is possible to thereafter obtain the parameter 
set P by the combinatorial optimization method. 

10 In the method of extracting physical model parameters of the third preferred 

embodiment, when it is judged that the number of updates of the parameter set P exceeds 
the predetermined number of times, the method of minimum-value solver is changed, 
instead of extracting the parameter set P at that time as one that gives the minimum value 
of the error function S, and verification on whether the function f y (v s , P) reproduces the 

15 observed values, judgment on whether the convergence condition for the minimum-value 
solver is satisfied and judgment on whether the parameter set P should be updated 
according to the number of updates are performed again. That makes it possible to 
adopt one of the combinatorial optimization method and the Newton-based solution as 
the minimum-value solver first and adopt the other one as the minimum-value solver to 

20 retry the extraction of the parameter set P when the first method fails the extraction of the 
parameter set P. 

While the invention has been shown and described in detail, the foregoing 
description is in all aspects illustrative and not restrictive. It is therefore understood that 
numerous modifications and variations can be devised without departing from the scope 
25 of the invention. 



